home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 November / macformat-018.iso / Demos / Extend 3.0 Demo / Demo Libraries / Demo Generic Lib / Demo Generic Lib.rsrc / MODL_1580_Input Random Number < prev    next >
Encoding:
Text File  |  1994-06-22  |  17.7 KB  |  824 lines
  1. real     p0[], p1[];        ** for plotter
  2. real     generalData[];
  3. real     saveMean,saveArg;    ** for dialog parameters when inputs used.
  4. integer    mySeed, tempSeed;
  5. integer numData;        ** found in solveGeneral()
  6. integer argK;
  7. integer oneInValid, twoInValid;                ** tests for input connectors being used
  8. integer    sending;
  9. real     values[], plotArray[], timeArray[];    ** for distribution plot
  10.  
  11. constant numMax is 50;    // size of data table (general)
  12.  
  13. **    This block generates random numbers
  14. **    Copyright © 1989-1992 by Imagine That, Inc.
  15. **    All Rights Reserved.
  16. **    Extend Generic Library, Input Random Number block; Alfy Riddle
  17. **        modified     2/5/92 JSL extensively modified for V2.0
  18. **                    9/1/92 BD modified Erlang
  19. **                    9/28/92 JSL added * 2.0 & factorB
  20. **                    9/30/92 JSL removed second get and setSeed from stepped
  21. **                    9/30/92 JSL reverted solveGeneral to 1.1 code.
  22. **                    2/16/93 JSL changed nummax to 50
  23. **                    5/4/93  JSL generate a number during initSim
  24. **                    10/14/93 PEH changed uniform real and intergers
  25. **                    1/17/94  JSL added tempSeed
  26. **                    1/26/94  DJK added triangular distribution
  27. **                    2/14/94  DJK modified trace and report
  28. **                    3/10/94  DJK added block label to report & trace
  29.  
  30. ** Distributions from Pritsker, Intro to Simulation & SLAM II
  31. ** Carroll, Simulation using Personal Computers
  32. ** and Gordon, System Simulation
  33.  
  34.  
  35. Procedure sortFcn() 
  36. {
  37.     integer i, numIn, x, y;
  38.     real    temp0, temp1;
  39.  
  40.         ** find the number of data points
  41.     if( noValue( data[0][0] ) )
  42.         return;
  43.     
  44.     i=0;
  45.     while( i<numMax && !noValue(data[i][0])  )    ** count Y values
  46.         i++;
  47.     numIn = i;
  48.     
  49.     for (x=0;x<numIn;x++)
  50.           for (y=x+1;y<numIn;y++)
  51.             if(data[y][0] < data[x][0])
  52.                 {
  53.                  temp0 = data[y][0];
  54.                  temp1 = data[y][1];
  55.                  data[y][0] = data[x][0];
  56.                  data[y][1] = data[x][1];
  57.                  data[x][0] = temp0;
  58.                  data[x][1] = temp1;
  59.                 }
  60. }
  61.  
  62.  
  63. **        these first three procedures are for the General distribution
  64. Integer getNumIn()
  65. {
  66.     integer i, numIn;
  67.     
  68.     i=0;
  69.     
  70.     ** find the number of data points
  71.     if( noValue( data[0][0] ) )
  72.         {
  73.         userError("No data to count in block number "+MyBlockNumber());
  74.         abort;
  75.         }
  76.  
  77.     while( i<numMax && !noValue(data[i][0])  )    ** count Y values 
  78.         i++;
  79.     numIn = i;
  80.  
  81.     if( stepped )    ** check for last two probabilities equal
  82.     {
  83.         **     IF the user is not interpolating and
  84.         **    the last two probabilities are not the same,
  85.         **    THEN we do not have a clear definition of a bin,
  86.         **    so we ASSUME a bin the same size as the last one
  87.         **    and add it to the list of points.
  88.         
  89.         if( realAbs( data[numIn-1][1]-data[numIn-2][1]) > 1.0e-15 )
  90.         {
  91.             data[numIn][1] = data[numIn-1][1];
  92.             data[numIn][0] = 2*data[numIn-1][0] - data[numIn-2][0];
  93.             numIn++;
  94.         }
  95.     }
  96.     
  97.     return( numIn);
  98. }
  99.  
  100.  
  101. Procedure solveGeneral()
  102. {
  103.     integer i;
  104.     real    factorA, factorB;
  105.     
  106.     **    this takes the density function (user input)
  107.     **    turns it into a cumulative distribution function
  108.     **    -> generalData[]
  109.     **    and scales the result to be a true CDF
  110.     
  111.     numData = getNumIn();
  112.  
  113.     **    would be best to sort the PDF ascending here
  114.     **    this will give fastest search later on
  115.     
  116.     **    "integrate" the PDF -> CDF
  117.     
  118.     makeArray(generalData, numMax);
  119.  
  120.     for (i = 0; i < numMax; i++)
  121.         generalData[i]    = 1.0;
  122.  
  123.     if( interpolate )        ** integrate trapezoidally
  124.         {
  125.         generalData[0] = 0.0;
  126.     
  127.         for( i=1; i< numData; i++ )
  128.             {
  129.             generalData[i] = generalData[i-1] + 
  130.                    0.5 * (data[i][1]+data[i-1][1]) * (data[i][0]-data[i-1][0]);
  131.             }
  132.         }
  133.     else if( stepped )    ** sum the probabilities * binWidths
  134.         {
  135.         generalData[0] = 0.0;
  136.     
  137.         for( i=1; i< numData; i++ )
  138.             {
  139.             generalData[i] = generalData[i-1] + data[i-1][1] * (data[i][0] - data[i-1][0]);
  140.             }
  141.         }
  142.     else    ** discrete
  143.         {
  144.         generalData[0] = data[0][1];
  145.     
  146.         for( i=1; i< numData; i++ )
  147.             generalData[i] = generalData[i-1] + data[i][1];
  148.         }
  149.  
  150.     for( i=0; i< numData; i++ )
  151.         {
  152.         data[i][1] *= 1.0 / generalData[numData-1];
  153.         generalData[i] /= generalData[numData-1];
  154.         }
  155. }
  156.  
  157.  
  158. Real findGen(real val)
  159. {
  160.     integer i;
  161.     real     result, slope, temp;
  162.     
  163.         **     first we find the bin of the data table (and CDF)
  164.         **        this random variable corresponds to.
  165.         **    Since the cumulative distribution function (CDF)
  166.         **        varies from 0 to 1 we can use a 0->1 real
  167.         **        random number to access our density function
  168.         **        VAL is assumed >= 0 and < 1
  169.  
  170.     i = 0;
  171.     slope = 0;
  172.     
  173.     while( val > generalData[i] )    ** find bin
  174.         {
  175.         i++;
  176.         }
  177.  
  178.     if( i > 0 )
  179.     {
  180.         if( interpolate )    ** since linear slope in PDF, analytic eqn found
  181.         {
  182.             slope = (Data[i][1]-Data[i-1][1]) / (data[i][0]-data[i-1][0]);
  183.                         
  184.             if (realAbs(slope) < 1e-10)    ** slope to zero, NaN result
  185.                 slope = 1e-10;            ** removes singularity at zero
  186.  
  187.             if( data[i][0] - data[i-1][0] == 0 ) ** interpolated fixed
  188.             {
  189.                 userError("Probability Range values in rows "+i+" and "+(i+1)+" are the same.");
  190.                 abort;
  191.             }
  192.             
  193.             temp = sqrt(data[i-1][1]^2 + 2.0*slope*(val-generalData[i-1]));
  194.             result = (-data[i-1][1] + slope*data[i-1][0] + temp)/slope;
  195.  
  196.         }
  197.         else if( stepped )        ** stepped
  198.         {
  199.             result = randomReal() * (data[i][0]-data[i-1][0])
  200.                         + data[i-1][0];
  201.         }
  202.         else
  203.             result = data[i][0];
  204.     }
  205.     else
  206.         result = data[0][0];
  207.     
  208.     return(result);
  209. }
  210.  
  211.  
  212. Real CalcValue()    ** for distribution plot
  213. {
  214.     real value, temp, min, max;
  215.     integer    i;
  216.  
  217.     tempSeed = RandomGetSeed();
  218.     RandomSetSeed(mySeed);
  219.     if( ranInt )        ** uniform integer
  220.         value = floor(randomReal()*(0.9999999999999+argDist-meanDist)+meanDist);
  221.     else if( ranReal )
  222.         value = (argDist-meanDist)*randomReal()+meanDist;
  223.     else if( norm )
  224.         value = gaussian(meanDist, argDist);
  225.     else if( binomial )    ** 11/28/90 BD Binomial bug workaround
  226.         value = RealAbs(DBinomial(meanDist, argDist));    ** take ABS value
  227.     else if( poisson )
  228.         value = DPoisson(meanDist);
  229.     else if( lognormal )
  230.         value = DLogNormal(meanDist, argDist);
  231.     else if( expo )
  232.         value = DExponential(1.0/meanDist);
  233.     else if( general )
  234.         value = findGen( randomReal() );
  235.     else if( erlang )
  236.     {
  237.         ** see Pritsker page714, note Pritsker's E{erlang}=U*K
  238.         
  239.         temp = 1.0;
  240.         for( i=0;i<argK;i++ )
  241.             {
  242.             value = RandomReal();    // from 0.0 to 0.99999....
  243.             
  244.             if (value == 0.0)        // if exactly zero, substitute 1e-3
  245.                 value = 1.0e-3;        // zero breaks erlang
  246.  
  247.             temp *= value;
  248.             }
  249.             
  250.         value = - meanDist * log(temp)/argDist;
  251.     }
  252.     else if( weibull )    // Pritsker p715
  253.         value = meanDist*(-log(randomReal()+1.e-6))^(1.0/argDist);
  254.     else if( hyper )
  255.     {
  256.         ** see Gordon, page 156, Note Gordon's E{hyper}=T
  257.         
  258.         if( randomReal() < argDist )    ** rate = 1/mean
  259.             value = DExponential(2.0*argDist/meanDist);
  260.         else
  261.             value = DExponential(2.0*(1.0-argDist)/meanDist);
  262.     }
  263.     else if (triag)                    // DJK 12/1/93 - triangular distribution
  264.     {
  265.         temp = randomReal();        // calc 0,1 random variable
  266.         min = meanDist;                // get the minimum
  267.         max = argDist;                // get the maximum
  268.         if(mode < min || mode > max)    // if the mode is out of range -> abort
  269.         {
  270.                 {
  271.                 Usererror("Most likely value is not between max and min for triangular distribution in Input Random Number block "+myBlockNumber()+".");
  272.                 abort;
  273.                 }
  274.         }
  275.             
  276.         if(temp <= ((mode - min)/(max-min))) // if on the left side of mode
  277.             value = min + Sqrt((mode-min)*(max-min)*temp);
  278.         if(temp > ((mode - min)/(max-min))) // if on the right side of mode
  279.             value = max - Sqrt((max-mode)*(max-min)*(1-temp));
  280.     }
  281.  
  282.     mySeed = randomGetSeed();
  283.     RandomSetSeed(tempSeed);
  284.  
  285. return(value);
  286. }
  287.  
  288.  
  289. procedure getNumber()
  290. {
  291.             **    let users have time varying statistics
  292.     if( oneInValid )
  293.         meanDist = oneIn;
  294.         
  295.     if( twoInValid )
  296.         argDist = twoIn;
  297.  
  298.     DistOut = CalcValue();    ** for distribution plot
  299.  
  300.     ** sysGlobal2 is the file reference number for the DEBUG TRACE
  301.     if( sysGlobal2 != 0.0 )    ** 0 is error, check for open file for TRACE
  302.         {
  303. // template for report:      |BLOCK NAME *****************|block number |BLOCK NUMBER*******
  304.         fileWrite(sysGlobal2,"Input Random Number          block number "+MyBlockNumber()+".  Current Time:"+currentTime+".","",True);
  305.         if(getBlockLabel(myBlockNumber()) != "")
  306.             fileWrite(sysGlobal2,"Block Label: "+getBlockLabel(myBlockNumber()),"",True);
  307.         if( oneInValid )
  308.             fileWrite(sysGlobal2,"     Parameter 1 = "+oneIn,"",True);
  309.         
  310.         if( twoInValid )
  311.             fileWrite(sysGlobal2,"     Parameter 2 = "+twoIn,"",True);
  312.             
  313.         fileWrite(sysGlobal2,"     Output = "+distOut,"",True);
  314.         fileWrite(sysGlobal2," ","",True);
  315.         }
  316. }
  317.  
  318.  
  319. on DistOut
  320. {
  321.     if (sending)
  322.         return;
  323.  
  324.     // if sysGlobalint9 == TRUE, msg came from plotter
  325.     // and if so, we should not recalculate the random number, 
  326.     // instead we should just let the connector value ride.
  327.  
  328.     if (!sysGlobalint9)
  329.         {
  330.         sending = TRUE;
  331.         if (oneInValid)
  332.             sendMsgToOutputs(oneIn);
  333.         if (twoInValid)
  334.             sendmsgToOutputs(twoIn);
  335.         sending = FALSE;
  336.         getNumber();
  337.         }
  338. }
  339.  
  340.  
  341. ** This message occurs for each step in the simulation.
  342. on Simulate
  343. {
  344.     getNumber();
  345. }
  346.  
  347.  
  348. integer    DataIsBad()        ** for distribution plot
  349. {
  350.     real temp;
  351.  
  352.     if( ranInt or ranReal )
  353.         if( meanDist > argDist )    ** swap if user input backwards
  354.         {
  355.             temp = meanDist;
  356.             meanDist = argDist;    ** the "a"
  357.             argDist = temp;        ** the "b"
  358.         }
  359.  
  360.     if( expo OR poisson OR lognormal)
  361.         if( meanDist < 1.0e-15 )
  362.         {
  363.             userError("Mean must be greater than zero in Input Random Number block number "+MyBlockNumber());
  364.             return(TRUE);
  365.         }
  366.  
  367.     if( binomial )
  368.         if( (meanDist < 1.0e-15) OR (meanDist > 1) )
  369.         {
  370.             userError("Probability must be between 0 and 1 in Input Random Number block number "+MyBlockNumber());
  371.             return(TRUE);
  372.         }
  373.             
  374.     if ( hyper AND (argDist < 0.0 OR argDist > 1.0))
  375.         {
  376.             userError("S must be between 0 and 1 in Input Random Number block number "+MyBlockNumber());
  377.             return(TRUE);
  378.         }
  379.     
  380. return(FALSE);    ** data is ok
  381. }
  382.  
  383.  
  384. ** If the dialog data is inconsistent for simulation, abort.
  385. on checkData
  386. {
  387.     oneInValid = oneIn;
  388.     twoInValid = twoIn;
  389.  
  390.     sysGlobal1 = 0.0;    ** prevent false reports
  391.     sysGlobal2 = 0.0;    ** prevent false debugs
  392.     
  393.     if (DataIsBad())    ** for distribution plot
  394.         abort;
  395. }
  396.  
  397.  
  398. procedure doInitSim()    ** for distribution plot
  399. {
  400.     sortFcn();
  401.  
  402.     if( general )
  403.         solveGeneral();
  404.     else if( erlang OR binomial )
  405.     {
  406.         argK = floor(argDist+0.5);    ** nearest integer
  407.         argDist = argK;
  408.     }
  409.  
  410.     saveMean = meanDist;    ** save dialog values in case inputs used
  411.     saveArg = argDist;        ** inputs write over dialog values
  412. }
  413.  
  414.  
  415. ** Initialize any simulation variables.
  416. on initSim
  417. {
  418.     if (randomSeed == 0)
  419.         mySeed = RandomGetSeed() + myBlockNumber()+1;
  420.     else
  421.         mySeed = randomSeed + myBlockNumber()+1;
  422.  
  423.     doInitSim();    ** for distribution plot
  424.  
  425.     if( getPassedArray(sysGlobal0, timeArray) )
  426.         getSimulateMsgs(FALSE);
  427.  
  428.     sending = FALSE;
  429.     getNumber();
  430. }
  431.  
  432.  
  433. string doName()
  434. {
  435.     string results;
  436.     
  437.         if( expo )
  438.             return "Exponential";
  439.         else if( erlang )
  440.             return "Erlang";
  441.         else if( weibull )
  442.             return "Weibull";
  443.         else if( ranInt )
  444.             return "Uniform Integer";
  445.         else if( ranReal )
  446.             return "Uniform Real";
  447.         else if( norm )
  448.             return "Normal";
  449.         else if( general )
  450.             return "General";
  451.         else if( poisson )
  452.             return "Poisson";
  453.         else if( binomial )
  454.             return "Binomial";
  455.         else if( lognormal )
  456.             return "LogNormal";
  457.         else if( hyper )
  458.             return "Hyperexponential";
  459.         else if( triag )
  460.             return "Triangular";
  461. }
  462.  
  463.  
  464. on endSim
  465. {
  466.     ** sysGlobal1 is the file reference number for the TEXT REPORT
  467.     if( sysGlobal1 != 0.0 )    ** 0 is error, check for open file for REPORT
  468.         {
  469. // template for report:      |BLOCK NAME *****************|block number |BLOCK NUMBER*******
  470.         fileWrite(sysGlobal1,"Input Random Number          block number "+MyBlockNumber(),"",True);
  471.         if(getBlockLabel(myBlockNumber()) != "")
  472.             fileWrite(sysGlobal1,"Block Label: "+getBlockLabel(myBlockNumber()),"",True);
  473.         fileWrite(sysGlobal1,"     Function = "+doName(),"",True);
  474.         if( general )
  475.             {
  476.             if( discrete )
  477.                 fileWrite(sysGlobal1,"       Discrete","",True);
  478.             else if( stepped )
  479.                 fileWrite(sysGlobal1,"       Stepped","",True);
  480.             else if( interpolate )
  481.                 fileWrite(sysGlobal1,"       Interpolated","",True);
  482.             }
  483.         fileWrite(sysGlobal1,"     Parameter 1 = "+meanDist,"",True);
  484.         if( oneInValid )
  485.             fileWrite(sysGlobal1,"     Parameter 1 input used","",True);
  486.         
  487.         fileWrite(sysGlobal1,"     Parameter 2 = "+argDist,"",True);
  488.  
  489.         if( twoInValid )
  490.             fileWrite(sysGlobal1,"     Parameter 2 input used","",True);
  491.  
  492.         if(triag)
  493.             fileWrite(sysGlobal1,"     Most Likely Value = "+mode,"",True);        
  494.  
  495.         if( comments != "" )
  496.             fileWrite(sysGlobal1,"     Comments = "+comments,"",True);        
  497.  
  498.         fileWrite(sysGlobal1," ","",True);
  499.         }
  500. }
  501.  
  502.  
  503. on createBlock
  504. {    
  505.     meanDist    = 0.0;
  506.     argDist        = 1.0;
  507.     expo        = 0;
  508.     erlang        = 0;
  509.     hyper        = 0;
  510.     weibull        = 0;
  511.     general        = 0;
  512.     ranInt        = 0;
  513.     ranReal        = 1;
  514.     norm        = 0;
  515.     poisson        = 0;
  516.     binomial    = 0;
  517.     lognormal    = 0;
  518.  
  519.     discrete    = 1;
  520.     interpolate = 0;
  521.     stepped        = 0;
  522.     
  523.     NMembers    = 200.0; ** for distribution plot
  524. }
  525.  
  526.  
  527. Procedure setText(string mReplace, string argReplace, string genReplace)
  528. {
  529.     mText = mReplace;
  530.     argText = argReplace;
  531.     genText    = genReplace;
  532. }
  533.  
  534.  
  535. on dialogOpen
  536. {
  537.     if( weibull )
  538.     {
  539.         setText("Scale =","Shape =","For General only");
  540.     }
  541.     else if( hyper )
  542.     {
  543.         setText("Mean =","s =","For General only");
  544.     }
  545.     else if( ranInt )
  546.     {
  547.         setText("Min =","Max =","For General only");
  548.     }
  549.     else if( ranReal or triag)
  550.     {
  551.         setText("Min =","Max =","For General only");
  552.     }
  553.     else if( general )
  554.     {
  555.         setText("unused.","(unused)","General values =");
  556.     }
  557.     else if( binomial )
  558.     {
  559.         setText("Prob =","N =","For General only");
  560.     }
  561.     else if( norm OR lognormal )
  562.     {
  563.         setText("Mean =","Std Dev =","For General only");
  564.     }
  565.     else if( expo OR poisson )
  566.     {
  567.         setText("Mean =","(unused)","For General only");
  568.     }
  569.     else if( erlang )
  570.     {
  571.         setText("Mean =","k =","For General only");
  572.     }
  573.     else
  574.     {
  575.         setText("Mean =","Arg =","For General only");
  576.     }
  577. }
  578.  
  579.  
  580. on weibull
  581. {
  582.         setText("Scale =","Shape =","For General only");
  583. }
  584.  
  585.  
  586. on hyper
  587. {
  588.         setText("Mean =","s =","For General only");
  589. }
  590.  
  591.  
  592. on erlang
  593. {
  594.         setText("Mean =","k =","For General only");
  595. }
  596.  
  597. on expo
  598. {
  599.         setText("Mean =","(unused)","For General only");
  600. }
  601.  
  602. on binomial
  603. {
  604.         setText("Prob =","N =","For General only");
  605. }
  606.  
  607. on poisson
  608. {
  609.         setText("Mean =","(unused)","For General only");
  610. }
  611.  
  612. on lognormal
  613. {
  614.         setText("Mean =","Std Dev =","For General only");
  615. }
  616.  
  617. on ranInt
  618. {
  619.         setText("Min =","Max =","For General only");
  620. }
  621.  
  622. on ranReal
  623. {
  624.         setText("Min =","Max =","For General only");
  625. }
  626.  
  627. on triag
  628. {
  629.         setText("Min =","Max =","For General only");
  630.         if(novalue(mode))
  631.             {
  632.             if(!novalue(meanDist) && !novalue(argDist))
  633.                 mode = (argDist - meanDist)/2 + meanDist;
  634.             }
  635. }
  636.  
  637.  
  638. on norm
  639. {
  640.         setText("Mean =","Std Dev =","For General only");
  641. }
  642.  
  643. on general
  644. {
  645.         setText("unused","(unused)","General values =");
  646. }
  647.  
  648. on discrete
  649. {
  650.     setText("(unused)","(unused)","General values=");
  651.     general = TRUE;
  652. }
  653.  
  654. on stepped
  655. {
  656.     setText("(unused)","(unused)","General values=");
  657.     general = TRUE;
  658. }
  659.  
  660. on interpolate
  661. {
  662.     setText("(unused)","(unused)","General values=");
  663.     general = TRUE;
  664. }
  665.  
  666.  
  667. on data
  668. {
  669.     if (!general)
  670.         userError("The General distribution must be selected to use this table.");
  671. }
  672.  
  673.  
  674. on doPlot
  675. {
  676.     integer i, numPlot;
  677.     real pStart,pEnd;
  678.     real pMax,pMin;
  679.     
  680.     sortFcn();
  681.  
  682.     solveGeneral();
  683.  
  684.         **    NOTE that the probability function may be interpolated
  685.         **        (i.e. the probability is linearly interpolated
  686.         **        between data points) or not. If it is
  687.         **        not interpolated the plotted points must define
  688.         **        each corner of the discrete bins (upper & lower
  689.         **        bound). This means twice as many points are
  690.         **        needed.
  691.         
  692.     if( interpolate )
  693.         numPlot = numData;
  694.     else if( stepped )
  695.         numPlot = 2 * numData - 1;
  696.     else    ** discrete
  697.         numPlot = 3 * numData;
  698.     
  699.     makeArray(p0, numPlot);        ** set up all arrays
  700.     makeArray(p1, numPlot);
  701.  
  702.     for( i=0;i<numData;i++ )        ** get values of function
  703.     {
  704.         if( interpolate)
  705.         {
  706.             p0[i] = data[i][0];    ** x coordinates
  707.             p1[i] = data[i][1];
  708.         }
  709.         else if( stepped )    ** use twice as many points so rect bins defined
  710.         {
  711.             p0[2*i] = data[i][0];    ** x coordinates
  712.             p1[2*i] = data[i][1];
  713.             
  714.             if( i < numData-1 )    ** odd number of points (skip last)
  715.             {
  716.                 p0[2*i+1] = data[i+1][0];    ** x coordinates
  717.                 p1[2*i+1] = data[i][1];
  718.             }
  719.         }
  720.         else    ** discrete
  721.         {
  722.             p0[3*i] = data[i][0];    ** x coordinates
  723.             p1[3*i] = 0.0;
  724.             p0[3*i+1] = data[i][0];    ** x coordinates
  725.             p1[3*i+1] = data[i][1];
  726.             p0[3*i+2] = data[i][0];    ** x coordinates
  727.             p1[3*i+2] = 0.0;
  728.         }
  729.     }
  730.  
  731.     pMax = 1;    ** the use of ceil & floor gives integer data limits
  732.     pMin = 0;
  733.     pStart = p0[0];
  734.     pEnd = p0[numPlot-1];
  735.  
  736.     installAxis(0, "", "Bin Range", FALSE, pStart,pEnd, "",
  737.         0, pMin, pMax, "", 0, 0.0, 0.0, blackpattern, blackcolor, numPlot);
  738.  
  739.     installArray(0, 0, "Range", p0, pStart,pEnd, 
  740.             numPlot, 0, blackPattern, redColor);
  741.     installArray(0, 1, "Density", p1, pStart,pEnd, 
  742.             numPlot, 0, blackPattern, redColor);
  743.  
  744.     makeScatter(0,0);
  745.  
  746.     showPlot(0, "Probability Density");
  747. }
  748.  
  749.  
  750. on plotNMembers
  751. {
  752.     real    min, max, val, points, range, pOverM;
  753.     integer    imembers, ipoints, i, j;
  754.     
  755.     if (randomSeed == 0)
  756.         mySeed = RandomGetSeed() + myBlockNumber();
  757.     else
  758.         mySeed = randomSeed + myBlockNumber();
  759.  
  760.     if (DataIsBad())
  761.         abort;
  762.         
  763.     if (novalue(NMembers) or NMembers <= 2.0)
  764.         {
  765.         usererror("Enter positive value for N greater than 2 (try 200)");
  766.         abort;
  767.         }
  768.     
  769.     doInitSim();    ** for distribution plot
  770.     
  771.     ** plot the distribution
  772.     makeArray(values, NMembers);
  773.     min = 1e4000;
  774.     max = -1e4000;
  775.     iMembers = NMembers;
  776.     for (i=0; i<iMembers; i++)
  777.         {
  778.         val = CalcValue();
  779.         
  780.         if (val < min)
  781.             min = val;
  782.         else if (val > max)
  783.             max = val;
  784.         values[i] = val;
  785.         }
  786.     
  787.     if (min == 0.0)
  788.         min = -0.01;
  789.     else
  790.         min = min-realabs(min*0.01);
  791.         
  792.     if (max == 0.0)
  793.         max = 0.01;
  794.     else
  795.         max = max+realabs(max*0.01);
  796.     
  797.     range = max-min;    ** correct for index
  798.     points = 50.0;
  799.     pOverM = 100.0/NMembers;
  800.     iPoints = points;
  801.     makeArray(plotArray, points);
  802.     for (i=0; i<iPoints; i++)
  803.         plotArray[i] = 0.0;
  804.     for (i=0; i<iMembers; i++)
  805.         {
  806.         j = (values[i]-min)*points/range-0.5;
  807.     
  808.         if (j < 0)
  809.             j = 0;
  810.         if (j >= ipoints)
  811.             j = ipoints-1;
  812.         plotArray[j] += pOverM;
  813.         }
  814.         
  815.     disposePlot(1);
  816.     installAxis(1, "Distribution Plotter", "Member Value", FALSE, min, max,
  817.             "Percent Members", FALSE, 0.0, 1.0, "", 0, 0.0, 0.0, 
  818.             blackpattern, blackcolor, ipoints);
  819.     installArray(1, 0, "% Members", plotArray, min, max, 
  820.                     ipoints, 0, blackPattern, cyanColor);
  821.     plotSignalFormat(1, 0, 1, 0);
  822.     autoscaley(1);
  823.     showPlot(1, "Distribution Plotter");
  824. }